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Abstract 



We use a previously introduced mapping between the continuum percola- 
tion model and the Potts fluid (a system of interacting s-states spins which 
are free to move in the continuum) to derive the low density expansion of 
the pair connectedness and the mean cluster size. We prove that given an 
adequate identification of functions, the result is equivalent to the density 
expansion derived from a completely different point of view by Coniglio et al. 
[J. Phys A 10, 1123 (1977)] to describe physical clustering in a gas. We then 
apply our expansion to a system of hypercubes with a hard core interaction. 
The calculated critical density is within approximately 5% of the results of 
simulations, and is thus much more precise than previous theoretical results 
which were based on integral equations. We suggest that this is because in- 
tegral equations smooth out overly the partition function (i.e., they describe 
predominantly its analytical part), while our method targets instead the part 
which describes the phase transition (i.e., the singular part). 
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I. INTRODUCTION 



Two previous papers in this series hereafter referred to as I and II, presented a gen- 
eral formalism of continuum percolation [0], where the system consists of classical particles 
interacting through a pair potential f(rj,rj), such that they can also bind (or connect) to 
each other with a probability p{fi,fj). Such a model is useful to describe microemulsions 
0, composite materials 0] or some properties of water |^. In this model, the clustering 
depends on the density p of the particles. As the density increases, so does the mean cluster 
size S. At a well defined critical density pc, S diverges, which signals the appearance of 
an infinite cluster. This is the percolation phase transition. Its critical behavior (i.e., the 
critical exponents) seem to be identical to the behavior of lattice percolation (see however 
a recent work by Okazaki et al. which claims differently) but the percolation threshold - 
the critical density - is sensitive to all the details of the system. 

The early theoretical attempts to calculate the percolation threshold culminated with 
the introduction by Balberg et al. of the notion of a critical total excluded volume Be 
101 . The excluded volume is the volume around one particle of the system in which the 
center of a second particle must be in order for the two particles to be connected. The 
total excluded volume of the system is therefore another measure of the number of particles 
in the system, i.e., of the density. When measured in this way, the numerical value of 
the percolation threshold seemed to be relatively insensitive to the shape of the particles 
(unlike the density itself). It was therefore considered an approximately universal quantity. 
However, this concept could not take into account other properties of the system, such as 
interactions between the particles, nor could it explain the remaining dependence of the 
percolation threshold on the details of the system. 

In 1977, Coniglio et al. proposed a different approach, based on a density expansion of 
the mean cluster size §]. This in turn served as the basis of an approximate calculation 
based on some integral equations analogous to those encountered in the theory of liquids . 
This approach yielded finally a theoretical prediction of the percolation threshold, but the 
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results remained only qualitative, discrepancies of up to 40 % with computer simulations 
being common [ p!0| , p!l[] . 

Recently, however, Alon et al. [|1^] and Drory et al. [jl3[ obtained quantitatively adequate 



results from Coniglio's expansion, by using the density expansion directly instead of integral 
equations. This posed a curious problem, because the critical density is not low enough to 
suggest that a power expansion should work. On the other hand, the expansion of Coniglio 
et al. had been developed to describe physical clustering in a gas, and its extension to 
general continuum percolation models was based on extensive analogies. It seemed that a 
more solid theoretical foundation was needed for continuum percolation before this puzzle 
could be addressed. Such a foundation has been laid in I and II, and we can now treat the 
problem of density expansions from a fresh point of view. The formalism presented in I and 
II is based on a quantitative mapping between the percolation model and an extension of 
the Potts model, the Potts fluid. For easy reference, we recall here the essential definitions 
and results. 

The s-state Potts fluid is a system of classical spins interacting with each other 

through a spin-dependent pair potential Vifi, Xf, Vj, Xj), such that 

{U (fi, r,) if Xi = A,- 
(1.1) 
W{f,,fj) if K^Xj. 

The spins are coupled to an external field h{r) through an interaction Hamiltonian 

N 

Hir,t = -Y.i;{X,)h{fi), (1.2) 

i=l 

where 

I s - 1 if A = 1 

m = (1-3) 

[ -1 if A 7^ 1. 

Up to some unimportant constants, the Potts fluid partition function (more precisely, 
the configuration integral) is 



^ = 1^ II / dfi---dfN exp 



i>j i=l 



:i.4) 



The magnetization of the Potts fluid is deflned as 



1 ^InZ ,^ , 



where h is the now constant external field. The susceptibility is 



dM , , 



N 

-/3^y(z,j)-/3EM^)V'(Ai 

i>j i=l 



The n-density functions of the Potts fluid are deflned as 

(fi, Ai; . . . ; fn, X„) = z(N^-n)\ J ^^"^^ ' ' ' d^N exp 

(1.7) 

Of particular interest is the spin pair-distribution function, g^'^\x,a;y,'f), deflned as 

g'^^\x,a;y,-f) = ^ p'^^\x,a;y,-f), (1.8) 
p{x)p{y) 

which tends to 1 when \x — y\ oo. Here p{x) is the local numerical density. It is often 
useful to deflne a spin correlation function h^'^^x, a;y,'y) as 

(f , a; y, 7) = gi'^ (f , a; y, 7) - 1. (1.9) 

This function tends to zero when \x — y\ — > 00. 

Any continuum percolation model deflned by v{i,j) and p{i,j) can be mapped onto an 
appropriate Potts fluid model with a pair-spin interaction deflned by 

U(i,j) ^v(i,j), 

exp[-(3W{t,j)] =q{t,j)exp[-(3v{t,j)], (1.10) 

where 

q{fi,rj) = 1 -p{fi,fj) (1.11) 

is the probability of disconnection . 

The relation between the Potts magnetization and the percolation probability, P{p), is 
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lim lim \im M = P(p). (1.12) 

/i— >0 Af— >oo 

For densities lower than the critical density pc, the susceptibility is directly related to 
the mean cluster size, S. 

lini lim Mm x = PS (p < Pc)- (1-13) 

ft— >0 N^cx) s^l 

An important quantity in the percolation model is the pair connectedness function 
g'^{x,y), the meaning of which is 

p{x) p{y) g\x , y) dxdy = Probability of finding two particles in regions 

dx and dy around the positions x and y, such 
that they both belong to the same cluster. (1-14) 

This function is related to the mean cluster size by 

S =l + p J dfg\f), (1.15) 

where we assume, as we usually shall, that the system is translationally invariant, so that 
g^(x, y) =g\x- y) and p{x) = p{y) = p. 

The pair connectedness is related to the Potts pair correlation functions by 

g\x, y) = lim (f, a; y, a) - g'f^ (f , a; y, r]) 



lim 



hf\x,a;y,a) - hf\x,a;y,r])] . (1.16) 



where the spins a and rj are arbitrary except for the conditions a, r/ 7^ 1 and cr ^ rj. 

Recently, Drory has applied this formalism to a non-trivial one dimensional model and 



managed to obtain its exact solution Continuing the investigation of the usefulness 

of this formalism, we consider in the present work the region p < pc and derive series 
expansion in powers of the density for the mean cluster size and the pair connectedness 
(the percolation probability, on the other hand, vanishes identically for these densities). 
To do this we introduce in section 11 the spin-functional differentiation, which is a simple 



generalization of the usual functional differentiation. With this tool, we obtain easily a 
diagrammatic expansion of the relevant quantities in section III. Section IV compares this 
expansion to the one derived by Coniglio et al. from a completely different starting point . 
Section V then applies the general results to a specific case, the extended hypercube models. 
The results are presented in section VI. Finally in section VII, we discuss the reasons for 
the quantitative success of the present approach, compared with the inadequacy of previous 
analytical attempts. 



II. SPIN-FUNCTIONAL DIFFERENTIATION 

The mathematical properties of the Potts fluid are very similar to the corresponding 
ones for a classical fluid. In the following derivations, we follow therefore very closely the 
presentation of Hansen and McDonald of the theory of classical fluids Occasionally 
we shall skip some mathematical details which are identical for the Potts fluid and for the 
classical one. 

Since all the quantities in which we are interested are expressible as statistical averages, 
we may choose to work in the grand canonical ensemble rather than in the canonical one. 
The grand canonical partition function, S, is given by 

OO 1 N N 

^=T.^ di-.-dN y: n^*(^'^«)n«^p[-w^'^')], (2.1) 

N=0 ■ '' Ai,...,Ajv j=l i>3 

where 

,.2\ 3/2 

m 



z*{i, \) = ( ) exp + mm^)] (2.2) 



is the generalized activity. Here, m is the mass of the particles (the spins), fi is the chemical 
potential, and h is Planck's constant. 

The n-density functions are now deflned to be 



^ 00 ^ N N 



p(")(fi,Ai;...;r„,A„) = ^ — — — y / rff„+i---df7v '[[z*{i,Xi)Y[exp[-(3V{i,j) 

^j^^Q[ly n). J A„+i,...,Aivi=l j>i 

(2.3) 



Paper II presented a generalization of the functional derivative, which will be called 
hereafter the spin-functional derivative. Let T\t{f, A)] be a functional of the function t(r, A), 
which depends on a position variable r as well as on an associated discrete spin variable A. 
Then the spin-functional derivative /5t is defined through the relation 

- j ir^j^^W.X), (2.4) 

where 5T is the change in T associated with a variation 5t in t(r, A). The only difference 
with the usual functional derivative is in the added summation over the spin variable. It is 
easily seen that this does not change any of the basic properties of the functional derivative 
operator. In particular, we have (see, e.g., Hansen and McDonald |^) that 

6t{x, a) 
St{y, 7) 

and 



S{x-y)5a,-y, (2.5) 



jdzY^t{z,\) 



Stiy, 7) 

The equivalent of the change of variable formula is now 



t(^,7). (2.6) 



6u{x, a) J 6v{x, A) 6u{x, a 

By a generalization of Eq. ( p.6|) , we now have that 



SE ^ N 



N N 



^ \ r 

.^.(^^. T.J^Jd2...dN E l[z*{^^,)l[exp[-m^,J)], (2-8) 

V ' / A''=l ■ A2,.--i-^jv *=2 i>j 

where fi = x and Ai = a. Therefore, comparing with the definition of the n-density 
functions, Eq. (p. 3D, we have that 

/ix,^ , z*(x,a) 5E ^.^ , 51nS . , 

^ ^ H bz*[x,a) ^ ' bz*[x,a) ^ ' 

An immediate generalization yields 



In particular, combining Eqs. ( p.9|) and ( |2.1CI| ), we have that 



51nH 



5z*{x, a)6z*{y, 7) 

A useful identity is obtained by using Eq. (p.5[) in conjunction with Eq. (p.9|), 



(2.11) 



Sp^^\x, a) 
6\n[z*{y,'y)] 



z*iy,7) 



6z*{y,'y) 



z*{x, a) 



6\nE 

6z*{x, a) 



p^^\x, a)5{x- y) 5^,^ + p^'^\x, a)p^^\y, -f)h^'^\x, a; 7), 



(2.12) 



where we have used Eqs. (|2.11|) , (|1.8| ) and (|1.9|) . 

By analogy with fluid systems , let us now define the spin direct correlation function, 

c(x,a;y,7), as 



c(f , a; y, 7) 



5 In 



p^^\x, a)/z*{x, a) 



M^H^>7) 

Then, with the identity Eq. ( p.5| ), we have that 

^ln[2:*(x,a)] 1 r - n 

5pW[y^^) p(i)(x,a) 

Finally, from the change of variable formula, Eq. ( p.7| ), we have that 

~5\n[z*{x,a)] _ r ^^~5\n[z*{x,a)] ~5p^^Kz,\) 



[2.13) 



(2.14) 



5{x-y)5a. 



^•^ 51n[z*(y,7)] 7 ^'^^ M'H^,A) 5 In [2;* (y, 7)]' 
Substituting Eqs. (|l|) and (|A|) into Eq. ( glSD , we find 



(2.15) 



h^^\x,a]y,-f) = c{x,a]y,-f) + / cii*^ p^^^(£. A) c(f, a; i*. A) /i^^)(f. A; y, 7), 



(2.16) 



which is the equivalent, for the Potts fluid, of the classical Ornstein-Zernike relation (see, 
e.g., Hansen and McDonald 0). 

To find out the imphcations of this for the percolation system, let us substitute Eq. ( p.l6| ) 
into the definition of the pair connectedness, Eq. ( |1.9| ). Then, we have that 

h^'^\x, a; y, a) — h^'^\x, a; y, 7) 

= c(x, a; y,a) + J dzp^^\z, a) c{x, a; z, a) h^'^\z, a; y, a) 
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+ / dzy^^ P^^\^, A) c(x, g; z, A) h^'^\z, A; y, a) — c{x, a; y, 7) 

— J dz p^''\z, a) c{x, a; z, a) h'^'^\z, a; y,^) — j dzp^^\z, 7) c(x, a; z, 7) h'''^\z, 7; 7) 

- 1 rf£^p«(^,A)c(f,«;£,A)/i(2)(£,A;y,7). (2.17) 



We use now the fact that for p < Pc the symmetry of the system is unbroken and therefore, 



p«(f,A) = -p(z), 
s 

where p{z) is the local density at z. For the same reason, we also have that 

c{x,a;y,a) = c(f,7;y,7), 
h^'^\x,a; y, a) = /i^^^ (x, 7; 7) , 
h^'^\x,a;y,-f) = h^'^\x, a;y, X), 



(2.18) 



(2.19) 



for any spins 7, A 7^ a, and other similar relations. Then we can rewrite Eq. ( 2.17|) as 



h^'^\x, a; y, a) — h^'^\x, a; y, 7) 

1 f 

= c{x, a; y, a) — c{x, a; 7) H — / dzp{z) c{x, a; z, a) h^'^^ {z, a; y, a) 



H / dz p{z) c{x, a; Zj'y) h^'^\z,'y;y,'y) / dzp{z)c{x,a;z,a)h'''^\z,a;y,'y) 

s J s J 

If s — 2 /■ 

dzp{z)c{x,a;z,'j)h^'^\z,a]y,a) / dz p{z) c{x,a] z,'-f) h^'^\z,'-f;y,'j). (2.20) 



Finally, we can take the limit s — > 1. Equation (p^.20D then becomes 



lim 

s^l 



h^^^ (x, a;y,a) - h^'^^ (f , a; y, 7) 
= lim I [c{x, a; y, a) - c{x, a; y, 7)] 

+ / dzp{z) c{x, a; z, a) — c{x, a; z, 7) h'''^\z, a; y, a) — h^'^\z, a; y, 7) 



We can now define a percolation direct connectedness function, c^{x,y), as 



(2.21) 



c\x, y) = lim [c(f , a; y, a) - c(x, a; y, 7)] . 

s— >1 



(2.22) 



Equation ( |2.21|) then becomes, with the help of the definition of the pair connectedness, 
Eq. (CT) , 



g^[x, y) = c\x, y) + dzp{z) c\x, z) g\z, y), (2.23) 



which is the percolation analog of the Ornstein-Zernike relation. This relation will prove 
useful in the next section. 

III. DIAGRAMMATIC EXPANSION 

The density expansion of the pair correlation of a Potts fiuid is obtained by following the 
steps leading to the corresponding density expansion for a classical fiuid, once some simple 
generalizations have been made. The following follows closely, therefore, the derivation 
presented in Hansen and McDonald 0. 

First, let us define the Potts equivalent of the Mayer /-function, as 

(f){fi,Xi;fj,Xj) = (t){i,i) = exp [-/3\/(r;, Ai;r;-, Aj)] - 1. (3.1) 
Let us now denote 

/(r^, Tj) = exp[-(3v{fi, fj)] - 1, (3.2) 

f*{fi, fj) = q{fi, fj) exp[-(3v{fi, fj)] - 1. (3.3) 



Then we have from the definition of V{i,j), Eqs. (1J-) and ( p,.l(]| ), that 



{ri,Xi;rj,Xj) = ■{ (3.4) 
f*{fi,fj) if Xi^Xj. 



The grand canonical function is now 

N N 



OO 1 „ iV iV 

5:— ycii-.-div y: E^*i^,>^^)E[^+<l^ihJ)]. (3.5) 



N=0 Ai,...,Aiv i=l i>j 

This expression is best represented diagrammatically, as 
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+ / . + A + A + - f^-^) 



In these diagrams, each circle corresponds to a function z*{i, Xi), and is therefore associated 
with a discrete spin variable as well as with a position coordinate. Each line corresponds to 
(f){i,j). For each black circle, we integrate over the space coordinate and sum over the spin 
coordinate. For example, 



• • 



J dfi dfj Yl ^i) ^j) 'Pi^i^ ^i'^ ^i' ^j)- (3-7) 



We will also have white circles over which there is neither integration nor summation. Note, 
finally, that the symmetry factor of the diagrams, being a purely combinatorial quantity, is 
the same here as for the usual diagrams without spin variables. While this is not necessarily 
so when the lines can stand for different functions as they can here, nevertheless it remains 
true because the identity of the lines is uniquely determined by spin variables which are then 
summed upon. 

It is now easy to see that all the usual definitions and theorems which hold for the 
usual spinless diagrams hold also for the spin diagrams introduced here. In particular, the 
Morita-Hiroike lemmas (see, e.g., Hansen and McDonald 0) hold here as well, provided the 
functional derivatives which appear in them are generalized to be the spin-functional deriva- 
tives introduced in the previous section. For example, we have the following generalized 
lemma. 

Lemma. If F is a diagram consisting of black t{f, A)-circles and (j) bonds, we have that 

-— — r = {all diagrams obtained by replacing a black t-circle of 
ot{x, a) 

F by a white 1-circle labeled (x, a)}. (3.8) 

From another lemma of Morita and Hiroike, we have immediately that (see Hansen and 
McDonald M), 



In S = {all simple connected diagrams consisting of black 2;*-circles and 0-bonds} 

= • + . . + /\ + ^A^ + (3.9) 

from which we can show that 

In p^^\x, a)/z*{x, a) = {all simple diagrams consisting of one white 

1-circle labeled {x, a), one or more black 
p(^)-circles and bonds, such that they are 
free of connecting circles}. (3.10) 

A connecting circle is a circle the removal of which (along with the bonds which emerge 
from it) causes the diagram to become disconnected. 

From Eq. (|2.13|) and the lemma ( p.8|) , we immediately obtain that 

c{x, a; y, 7) = {all simple diagrams that consist of two white 

1-circles labeled (x, a), and (y, 7), black p*^^^-circles 
and (f) bonds, and are free of connecting circles} 



+ 



(3.11) 



In all these diagrams, one of the white circles is labeled (x, a) and the other is labeled {y, 7). 

The Ornstein-Zernike relation, Eq. ( p.l6|) , now determines h^'^\ For vanishing magnetic 
field h = 0, we have significant simplifications, since p^^\z, X) = (l/s)p, where p = (N) /V 
is the usual density. Furthermore, h^'^\x, a;y,'j) and c{x,a;y,'j) now depend on x — ^ 
provided v{x, y) = v{x — y) and p(x, y) = p{x — y), which we always assume. In this case, 
the Ornstein-Zernike relation is very simple in Fourier space, i.e., 

h^'\k,a,,)- , f^-^;i> (3.12) 

1 - (p/s)c(fc,a,7) 

where 
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c{k, a, 7) = / rff e"^^(^'^^) c(x, a; y, 7). (3.13) 



Going now to the percolation picture, we notice that Eq. ( |1.15| ) for S can be rewritten 



S = l + pg^O). (3.14) 



From the definition of g\ Eq. ( |1.16| ), and from Eqs. ( p.23| ), ( p.22| ) and ( p.l2| ), we have now 



^ l-pct(O) 1 -plim [c(0,a,a) -c(0,a,7)]' ^^'^^"^ 

This equation, together with the expansion of c(x, a; y, 7), Eq. ( |3.11|) , defines an expan- 
sion of 5" in powers of the density. The critical density pc is the radius of convergence of this 
series, which we can find by using standard methods of convergence analysis (see, e.g., ref. 



IV. EQUIVALENCE WITH CONIGLIO DEANGELIS FORLANI EXPANSION 

In 1977, Coniglio, DeAngelis and Forlani (CDF) obtained a density expansion for 
c'^lxj'y) from a completely different point of view. We discuss now the connection of the 
present work to the CDF expansion. 

Coniglio et al. start by considering a classical fluid (no spins), where the inter-particle 
interaction is v{i,j). For this fluid, the classical Ornstein-Zernike direct correlation function, 
Cf{x,y), has a density expansion given by Eq. ( 3.11| ), except that the white 1-circles are 



labeled just x and y and the bonds represent the Mayer /-function deflned as f{i,j) = 
exp[—(3v{i,j)] — 1, in the same notation used here in Eq. ( |3.3| ). Coniglio et al. then assume 
that the Mayer /-function can be decomposed into a sum f{i,j) = f\h 3) + f*{i, j) where p 
is related to the particles being connected and /* to their being unconnected. The speciflc 
definitions of /^ and /* used by Coniglio et al. are relevant only to physical clustering 
in a gas, but the assumption of the decomposition of f{i,j) was extended by analogy to 
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other systems by several workers P, [TO| , pl| . However, no formal basis was provided for the 



definition of and /* in the general case. If we take /* as formally defined here in Eq. ( p.3[ ), 
which is natural since this function is directly related to the probability of the two particles 
being unconnected, we find that 

f\^,J)^fi^,J)-n^,J) = pit, j) exp[-Pv{t, j)] . (4.1) 

This result is satisfying since turns out to be directly proportional to the probability of 
the two particles being connected. 

Coniglio et al. then replace / in the diagrammatic expansion of Cf{x,y), with the sum 
f'^ + f*, and consider all the decomposed diagrams thus obtained, in which every bond is 
either an /''" or an /* function. The sum of all the diagrams which contain at least one 
continuous path of functions between the white 1-circles, they consider to be c^{x,y), 
while the remainder is called c*(x, y), so that c/(x, y) = c^{x, y) + c*(x, y). 

It is hardly obvious now that this algorithm agrees with the expansion obtained here. 
Indeed, the function /''" never appears in the present expansion. Nonetheless, the two for- 
mulations are equivalent. 

The proof of the equivalence holds diagram by diagram. Therefore we can take for 

• • 

simplicity a specific example, e.g., the diagram : : . In this diagram , every dotted 

o o 

line denotes a function f{i,j). Then, if we consider all the diagrams obtained by replacing 
every /-bond by either an /*-bond or an /^-bond in all possible ways, and keep only those 
diagrams which contain a continuous path between the two white circles, we obtain 



CDF 



• • 

6 6 



• • 

(4.2) 



where a double line denotes a function /^(i,j), and CDF[diagram] on the left hand side 
means the result of the CDF algorithm when applied to the given diagram. The dotted 
lines on the right hand side (/-bonds) arise because once a continuous /^ path exists, all the 
diagrams obtained from further replacements of /-bonds by /^-bonds and /*-bonds must be 
counted. This sum can be expressed more compactly by keeping the original /-bonds. 
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We now wish to prove the equivalence of this scheme to the one we have developed here, 



i.e, 





• • 




CDF 




= lim 




6 6 





c — o 

a a 



o — o 
a Tj 



(4.3) 



where on the right hand side, every diagram is of the spin type, and every full line represents 
a function (b(i,j). To simplify still further, let us ignore the integration over the coordinates 
of the black circles (but not the summation over spins in the spin diagrams) . The integrals 
can then be added at the end of the proof. E.g., in this section we shall assume that 



• • 

6 6 



3 4 
• • 

6 6 

1 2 



/(l,3)/(3,4)/(4,2)/(l,2). 



(4.4) 



and 



o — o 



= ^ 0(1,3)0(3,4)0(4,2)0(1,2), (4.5) 
l,a 2,77 

where the labels 1, 2, 3, 4 denote the spatial coordinates of the vertices. Note that in the last 
diagram, Ai = a, A2 = are fixed, but we sum over the spins A3 and A4. 

With this convention, we are ready to prove the equivalence of the two formulations. To 



do this, let us define a new "bond percolation" problem on the four-vertex graph ^ ^ 

In this new problem, every bond can be either "open" (in which case its two ends are 
said to be connected) or closed (in which case the ends are disconnected). Let us further 
assume that the bond connecting the vertices i and j is open with a probability P{i,j) 
and closed with a probabihty f*{i,j). Formally these are not true probabilities since their 
sum is not normalized to 1, but this is of no importance. One can always normalize the 
probabilities by dividing by the proper product of / functions. With this convention, every 
diagram obtained from the CDF algorithm now represents the (unnormalized) probability of 
a specific configuration of open and closed bonds on the four-vertex graph. E.g, consider 



0***0 



/t(l,3)/t(3,4)/t(4,2)r(l,2), 
15 



(4.6) 



where a line of * represents a function /*. The same diagram can also be thought of as a 
configuration of open and closed bonds on the four-vertex graph, by letting a double line 
represent an open bond and a *-line a closed bond. This geometrical configuration occurs 
with a probability /''"(I, 3)/^(3, 4)/^(4, 2)/*(l, 2), which is just the value of the equivalent 
diagram, Eq. ( [4. 61 ) . Every diagram therefore functions in a double capacity. On the one hand 
it represents an actual configuration of open and closed bonds on some graph. On the other 
hand, it represents some functional value. This functional value is just the (unnormalized) 
probability of the bond configuration represented by the same graph. 

Now, two vertices which are connected by a path of open bonds can be said to belong 
to the same cluster. In particular, the (unnormalized) probability that the two white circles 
belong to the same cluster, denoted Pci{l,2) is the sum of all the diagrams in which a 
continuous path of open bonds exists between these circles. This however, is exactly the 
class of diagrams selected by the CDF algorithm. Hence, 



PJh2)=CDF 



Therefore, what we now need to prove is that 



• • 

6 6 



(4.7) 



P./ (1,2) 



lim 



o — o 



o — o 



a a 



(4.8) 



a rj 

The prove this we show, by paralleling the derivation given in I, that there exists a mapping 
of the new bond percolation problem onto an s-spin model. Then Pci{l,2) ( which plays 
here the part of the pair connectedness) is equal to the difference of two spin diagrams. The 
formal proof of this claim is derived in the Appendix. It is easily seen to hold for every 
diagram in the CDF expansion, and shows therefore that this expansion is indeed equivalent 
to the one derived here in the context of the Potts fluid. 
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V. APPLICATION TO EXTENDED HYPERCUBES 



Let us see now how the expansion obtained in section III can be apphed to finding the 
critical density pc- The system we consider consists of D dimensional hypercubes which 
have an impenetrable (hard) core of side a, surrounded with a hypercubic permeable shell 
of side (i > a, so that two particles are bound if their shells overlap. In other words, if 



' I ~ K-^i 1 ■ ■ ■ 1 -^t 



f ) and = (x], . . . , x'j) are the positions of two particles, then 

00 if for all 1< k < D, \x'f — x'^A < a 
v{i,3) = { - - ^' (5.1) 

otherwise. 

1 if for all 1 < k < D, \xf — x'^l < d 
P{^,J) = { " " ' (5.2) 



otherwise. 



Hence, from Eq. 



, -1 if for all 1 < k < D, \xf — x''\ < a 
fir^,r,) = { - - ^' (5.3) 

otherwise. 

-1 if for all 1 < k < D, Ix'f — x'^l < d 
nr^,f,) = { - - , I . ,1 ^^^^^ 

otherwise. 

We choose this rather unrealistic system because for any dimension D the required calcu- 
lations always turn out to be much simpler for cubes than for the more realistic spheres. 
Since the aim at present is to prove the usefulness of the method rather than investigate a 
specific system, one should not be overly troubled by this choice. Let us note, however, that 
the spherical version of this model has been used to model microemulsions 0, so that it is 
at least mildly relevant to some real systems. It is also the simplest model which contains 
interactions, and therefore a good testing case. 

Let us now define two functions of a single variable. 



1 if \x\ < a 

(5.5) 

if Ixl > a. 
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1 if Ixl < d 

F2{x) = { (5.6) 
if \x\ > d. 

The equations (|5.4|) can now be rewritten as 

f{x\...,x'') = -f[F,{x^). (5.7) 

i=l 

nx\...,x'') = -l[F2ix^). (5.8) 



i=l 



We shall now calculate up to the second order in p (i.e., up to square diagrams). The 
zeroth order is 



lim 



a a a 7 

where we have defined the aspect ratio, r], as 



I dz [fiz) - riz)] = (l - r^^) , (5.9) 



a 

r]^-. (5.10) 



For the first order (triangular diagrams), we need, e.g., the integral 



n{y)^ Jdxfix)f{x-y) = f[I^Jdx'F,{x')F,{x'-f)Y (5.11) 

It is precisely the factorization of this and similar integrals into a product of identical one 
dimensional integrals that is the main simplification afforded by the use of a cubic geometry. 
As a beneficial side-effect, it allows us to check the influence of system dimensionality on 
the critical density. 

fl{y) is very simple to evaluate when we realize that / dx Fi{x) Fi{x — y) is merely the 
overlap of two segments of length a, one of which is centered at the origin, the other centered 
at the position y. By direct inspection, we have therefore that 



r 2a — Ivl if \y\ < 2a 

J dxF,{x)F,{x-y) = l , (5.12) 



2a - 


\y\ 


if 


\y\ 


< 2a 







if 


\y\ 


> 2a 


2d- 


\y\ 


if 


\y\ 


< 2d 







if 


\y\ 


> 2d 



r — o 11 \y\ <^ z,a 

j dxF2{x)F2{x-y) = I , (5.13) 



J dx Fi{x) F2{x — y) 



2a if \y\ < d — a 
d + a — \y\ if d — a<\y\<d + a- (5-14) 
if \y\ > d + a 

Let us consider a typical contribution from the triangular diagrams, e.g., the integral 
/ dx dy F2{x) F2{x — y) Fi{y). In this integral, we interpret the term Fi{y) as merely spec- 
ifying the limits of integration, i.e., determining that \y\ < a. Since d > a, we have from 
Eq. (|1D that 



J dxdyF2{x) F2{x - y) Fi{y) = J dy J dx F2{x) F2{x - y) 



Ada — a^. 



(5.15) 



All the required one dimensional integrals can now be calculated along the same lines. In 
fact, equations ( |5.12| )-( |5.14| ) suffice to calculate all the triangular and square diagrams, 
except for the last fully connected square diagram which is a product of six F functions. 
For this diagram the previous arguments have to be generalized but it can be done easily 
along the lines already presented. 

All the integrals we need are presented in Tables I-IV. The first column in each table 
contains the notation of the integral, and the third column its value. All these integrals are 
Z^-tli powers of some one dimensional integrals, [see, e.g., Eq. ( ^.11|) ]. The corresponding 
one-dimensional integrals appear in the second column. 

We now have, after some straightforward calculations, that the contribution of the tri- 



angular diagrams, denoted {2df'^ k2, is, in the notation of Tables I-IV, 



{2df'' k2 = h-2l2 + h 



(5.16) 



Similarly, the contribution from the square diagrams is denoted {2df'^ k^^ and in the 
notation of Tables I-IV is equal to 



(2d)'^ k^ = ^ Ji - 5J2 + 5 J3 - - m + ^K2 - ^K, + 7K^ - lOi^s + 

+ ^Li-L2-L3 + ^L4-L5. (5.17) 

Substituting these results into Eq. (|3.15| ) yields, for the mean cluster size. 
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The critical density is commonly measured in dimensionless units. Let us define 

B = {2d)^p, (5.19) 

which is chosen to reduce to the total excluded volume when a ^ 0. As mentioned in 
the introduction, this reduces the dependence of the density on the details of the system 



and facilitates the presentation of the results. Finally, using standard methods [|16|, we find 
the series expansion for the mean cluster size to be 

S=l + SiB + S2B^ + S3B^ + 0{B^), (5.20) 

with 

52 = k2 + Si 

53 = k3 + 2Sik2 + (5.21) 



VI. RESULTS FOR HYPERCUBES 

To find the critical density, we need some extrapolation method which will yield the point 
of divergence of the series (|5.20|) . Unfortunately, we know very few terms in this series. The 
results can be improved, however, by using additional information. In particular, S diverges 
as S* ~ {Be — B)~^ , and 7 appears to be universal, i.e., independent of the interaction, and 
furthermore, equal to the well known value it takes in lattice percolation. Biased methods 
use the known value of 7 to calculate Be. There are several methods available, but for the 
present use, the best seems to be the one developed by Arteca, Fernandez, and Castro (AFC) 
|17| . The biased version of this method relies on the fact that the function {B^ — ByS{B) 
is analytical around Be- Therefore, we define a function of two variables. 
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giu,B)=[l-^y S{B), (6.1) 

where we have assumed a known value for 7. Given an expansion S = SnB^, we can 
expand g{u, B) in powers of 5, with the result, 



g{u,B) = Y.9n{u)B\ (6.2) 



n=0 

where 



^?n(«) = E(^^j(-«)"'^n-fc- (6-3) 

Since g{u,B) is analytical at u = Be, the idea is now to look for values of u which maximize 
the convergence of the series {gn{u)}'^^i as n — > 00. The AFC proposal is to generate a 
series {un}'^=i of solutions of the equation 

gn{Un) = 0. (6.4) 

Then we expect that 

lim Un = Be- (6.5) 



n— »oo 



In the case of the series ( [5. 201 ), we know the terms only up to n = 3, which is too little to 



extrapolate a reliable limit to the series {un}'^=i- Therefore, we take simply the last term of 
the series , M3, as our estimate of Be- 

In figures (l)-(5), we compare this theoretical estimate to results of computer simula- 
tions for dimensions 2, 3, 4 and 5, as a function of the aspect ratio rj. The simulations were 
performed with an efficient new algorithm we introduced recently for investigating contin- 



uum percolation with interactions [|T8[. Unlike the usual Metropolis algorithm, this method 
increases serially the density in the system until percolation is achieved. The effect of the 
interactions is included through a rejection criterion which produces an effective statistical 
weight for the configurations. The final (percolating) configuration contained close to 30,000 
particles in two and three dimensions and around 10,000 in the higher dimensions. The re- 
sults are averages over 10 independent runs and the numerical error is estimated to range 
from 5% to 10 %. 
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Let us first consider the well reproduced qualitative behavior of Bc{ri). The main features 
can be understood from simple arguments. As pointed out by Bug et al. 0, the minimum 
of Bc{ri) results from the competition between two processes. As the diameter of the hard 
core increases, it becomes harder to bring particles close enough to each other for them to 
bind. This effect clearly dominates at high rj. At low 77, on the other hand, the hard core 
is much too small relative to the soft shell to prevent binding significantly, but it increases 
the average distance between bound particles. As a result the clusters are "longer", and 
therefore percolate "sooner", hence Be is lower. 

The dependence on D is understandable in terms of the ratio of the permeable shell vol- 
ume to the hard core volume. Clearly, this ratio increases with the dimensionality. Therefore, 
at any given 77, the influence of the hard core diminishes as D increases. As a result, the 
graph of Bc{ri) looks flatter, while the influence of the final singularity at 77 = 1 (at which 
there is no possibility of binding anymore and therefore no percolating phase), is limited to 
regions of higher 7]. This influences also the feasibility of the simulations. Clearly, the larger 
the hard core, i.e., the greater 77, the harder it is to simulate the system. Hence, e.g., in 
two dimensions, we have no results beyond rj = 0.8 due to this difficulty. As D increases, 
however, higher rj become more accessible. Unfortunately, it is generally harder to simulate 
high dimensional systems, which is why we stopped at D = 5. The trends, however, are 
obvious, and there should be nothing exceptional for D > 5. 

Turning now to the detailed comparison between theory and simulations, we see that 
the results are quantitatively close to each other. Even in the worse case, in two dimensions 
(Fig. 1), the difference between theory and simulations is at most 11% (at t] = 0.5). The 
position of the minimum, theoretically predicted to be t] = 0.79, also agrees well with the 
simulations. 

Very recently, however, Okazaki et al. have performed new simulations of two dimen- 
sional continuum percolation, and claim that the critical exponents are in fact different from 
those of lattice percolation. For the exponent 7 they find 7 = 1.94 instead of the higher 
43/18 obtained in lattice percolation [jl9|. More work is certainly needed to resolve this 
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issue. Here we merely note the effects of such a shift in the value of 7 on the theoretical 
prediction. Figure 2 presents the results obtained in two dimensions when the value of 7 
found by Okazaki et al. is used for the biasing. The fit is quite improved by this change. Un- 
fortunately the approximations we use are too drastic to claim that our calculations actually 
support the value obtained by Okazaki et al.. But they do suggest that the question of the 
universality class of continuum percolation is not yet settled and requires more investigation. 

Note however that the relative change in Be is smaller than the change in 7 which 
produced it. Therefore, Be is not extremely sensitive to the value of the bias. Since other 
simulations find values of 7 in higher dimensions which are at least close to the values 
obtained in lattice percolation, we may use these lattice values for the bias and expect that 
Be will not be affected notably. For this reason and for lack of other data on the value of 7 
in continuum percolation, we therefore use the accepted lattice values for dimensions higher 
than two. 

In three dimensions, up to r] — 0.8, the greatest discrepancy between theory and "ex- 
periment" is only 4 %, well within the limits of the simulation error. Even at 77 = 0.8, the 
largest discrepancy we have, the deviation is only about 6 %, still within the simulation error 
(which increases slightly with 77 because of the above mentioned problems in performing the 
simulations in the vicinity of 77 = 1). 

This trend continues in dimension 4, where the discrepancy between the theory and the 
simulations remains below 6 % until 77 = 0.65, and in dimension 5. Note, however, that 
as a general rule simulation errors increase with dimensionality due to the greater difficulty 
in performing the simulations, and therefore comparison with the theory becomes slightly 
more problematic. In particular, Ski D — b there seems to be a systematic discrepancy with 
the theory which could well reflect some finite size scaling effects in the simulations. It 
is therefore quite possible that the theory is actually more precise at this point than the 
simulations. 

Considering that we use the series for S only up to second order, the quantitative agree- 
ment obtained is a remarkable success for the theory. 
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VII. DISCUSSION 



We have seen how the Potts fluid mapping allows us to derive a general expansion in 
powers of the density for the mean cluster size and the pair connectedness. This expansion 
can be used to calculate the percolation threshold by calculating the first terms in the series 
then using some extrapolation method to find the radius of convergence of the series. We 
have applied this scheme to interacting hypercubes in dimensions ranging from 2 to 5, and 
we have shown that even with merely three terms in the series we already obtain very good 
quantitative results, within a few percents of the results of the computer simulations. 

Ours is not the first work to use density expansions of the mean cluster size, obtained in 
one way or another (e.g., from an analogy with lattice systems [^), but it is the first one to 
derive quantitatively adequate results for interacting systems. We believe this success stems 
from a simple but fundamental reason. 

When Coniglio et al. [|| obtained their expansion of S', they did not use it directly to 
calculate the critical density. Instead, they turned for inspiration to the theory of liquids and 
considered the integral equations for the pair correlation which had proved successful there, 
primarily the Percus- Yevick (PY) equation. Coniglio et al. derived a percolative analog of 
this equation by applying the CDF algorithm to its diagrammatic expansion. Some time 
later, De Simone et al. |]TD[ applied this equation to a system of extended spheres, the spher- 
ical analog of the system investigated in the present work. The results were qualitatively 
correct but quantitatively poor. Discrepancies between theory and simulations ranged from 
40 % to no less than 10 %. Several researchers fill] extended this work to other interactions 
or shapes with always the same quantitatively inadequate results. We believe this failure 
holds a simple (in hindsight, in fact, obvious), but apparently under-appreciated lesson. 

The guiding principle behind the work of Coniglio et al. and their followers seems to have 
been viewing continuum percolation as a high density phenomenon. This is very reasonable 
from a point of view which starts from the usual theory of fiuids and then extracts the 
percolative quantities from their normal fiuid analogs. The CDF algorithm is explicitly 
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based on this point of view; given, e.g., the direct correlation function c/(r) of the fluid, the 
percolating part, c^(r) can be extracted from it by looking at a diagrammatic expansion. 
This means we need to know c/(r) first. The percolation transition occurs at relatively high 
densities, so that the normal fiuid is either a dense gas or a liquid. In this case, integral 
equations are the main tool for calculating the fiuid's properties, and therefore an adequate 
starting point for obtaining the percolative properties as well. 

But such a point of view misses the most important aspect of the percolation transition, 
namely that it is a phase transition. Integral equations are known to describe the liquid well 
only away from the liquid-gas critical point. As we approach the transition, they fail. The 
quantitative failure of the percolative integral equations suggests that they suffer from an 
analogous defect. There is direct evidence for this claim. Seaton and Glandt p2| checked 



the mean cluster size for a wide range of densities and compared the predictions of integral 
equations with the simulations. They found that the integral equations are indeed very 
successful away from the critical density, but worsen steadily as p pc- In short, the 
percolative integral equations behave with respect to the percolation transition exactly as 
their analog in liquids behave with respect to the liquid transition. 

At first sight this seems puzzling. At pc, the normal liquid is usually away from its own 
critical point, and it is therefore well described by integral equations. Now the procedure 
for extracting the percolative part from these equations involves no further approximations. 
Why then, do we obtain a relatively bad approximation by extracting exactly the percolative 
part from a very good approximation? 

This is precisely the heart of the matter, which is the trivial observation that the phase 
transition is entirely ruled by the singular part of the relevant functions. Integral equa- 
tions tend to smooth out this singular part, as evidenced, for example, by their systematic 
overestimation of the critical density [^. This doesn't matter for the description of the 
liquid, because the singular part is small at these densities. However when we turn to the 
percolative analog, the entire physics we are after lies precisely in this overly smoothed out 
singular part. 
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To understand a singular behavior, we must use methods adapted to singular functions. 
This is the simple but under-appreciated lesson from all the preceding. It requires us to 
abandon the view of continuum percolation as primarily a high density process. The density 
at which the transition takes place is irrelevant for the choice of descriptive tools, since it only 
influences the analytical part of the function. The proper tools to describe the percolation 
transition are not those which describe well high density systems, but rather those which 
describe well other phase transitions. 

The present work proceeds clearly from this point of view. The underlying normal liquid 
plays no part in the description presented here. Nowhere does the function c/(r) or any of 
its parents appear. Instead, we have a mapping onto a Potts fluid, from which we calculate 
all the relevant quantities. One important difference between this and the older point of 
view is that the magnetic transition of the Potts fluid can be identified with the percolation 
transition, since the two relevant order parameter are essentially the same. Therefore we are 
forced from the outset to employ only methods which are useful for the description of the 
magnetic phase transition. Indeed, integral equations are not a natural tool to use here since 
they will clearly fail in the description of the Potts fluid. Instead, we rely on one of the best 
tried methods in critical phenomena, namely, an expansion valid in the disordered phase 
(high temperature usually, low density in the present case). A low density expansion makes 
no sense from the point of view of percolation as a high-density phenomenon. Indeed, we 
do not expect the actual value of S obtained from the truncated series to be even remotely 
good. But what we look for is a failure of convergence, the appearance of a singularity, and 
for this a truncated series can be singularly well adapted, as countless investigations of the 
Ising model have shown repeatedly ||T5| . 

It is precisely because the mapping onto a Potts fluid identifies one phase transition 
(percolation) with another (magnetic) that we can achieve a remarkable quantitative agree- 
ment between theory and simulations. Clearly, the truncated series expansion coupled with 
an extrapolation method to find the radius of convergence is at present the best available 
method to calculate theoretically the critical density in continuum percolation. 
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APPENDIX 



• — 

We consider some diagram in the expansion, e.g., the example used in section 0, ^ 

On this diagram we define a bond-percolation problem. We can now map the problem onto 
a spin model by simply assigning to every vertex an s-spin, and assigning to every bond 
between the vertices i and j a function defined as in Eq. (|3.1|) . We now proceed 

to show, just as in paper I, that the bond-percolation model corresponds quantitatively 
to the spin model. The result is again a geometrical mapping according to which every 
percolation configuration corresponds to several spin configurations, in each of which all the 
spins belonging to a single cluster are parallel. Different clusters, however, have randomly 
assigned spin values. The following derivation parallels closely the one presented in paper I 
for the continuum case. 

Let the chosen diagram contain n vertices, and let us denote 

r(l,2,...,n)= n '^(^'i')- (Al) 

all existing 

bonds in 
the graph 

By convention, the two white circles are always labeled 1 and 2. The value associated with 
the diagram in the spin model is therefore I]a3,A4,...,a„ r(l, ...,«), where we assume that the 
spins of the two white circles, Ai and A2, are fixed , and where we ignore the integration as 
in section |V[ The integration can be restored at the very end of the proof without changing 
anything. For example, for the four-vertex graph mentioned above, we have 

r(l, 2, 3, 4) = 0(1, 3)0(3, 4)0(4, 2)0(2, 1). (A2) 

Let us now select two spins, say Aj and Aj such that the bond 4>{i,j) exists in the graph. 
In the expression for r(l, . . . ,n), let us separate all possible configurations {A^} into those 
where Aj = Aj and the rest. Then, 

27 



We can rewrite the sum over the spins when Aj 7^ as the difference between the sum over 
the spins without constraints and the sum over the spins when Aj — Xj, so that 

^r(l....,n)=fHiJ) E %7r + /'('.i)5:^^^i-f^. (A4) 

where we used the relation f^{i,j) — f{i,j) — f*{hj)- The last sum J2{x^} on the right 
hand side is now performed over all spin configurations without constraints. Let us now 
choose another pair of spins, say Aj and (if such a bond exists). Repeating the previous 
procedure, we obtain that 

This argument can be easily generalized to all pairs of spins. Let us consider one of the 
sums into which the function J2x^ r(l, . . . , n) has been decomposed a step before, and con- 
sider a pair (m, n) such that the bond 0(m, n) exists in the chosen diagram. Two possibilities 
arise: 

1. Previous constraints already determine that A^ = A„ (for example, there could be 
some p for which A^ = Xp and A„ = Ap). Then, 



Ari 



{ij)---(j){l,k) ( ^ . (l){i,j)---(l){l,k)(l){m,n) 



< previous > < prcv. > 

I constraints J I constr. J 

+ r(m,n) V .,^^^v;/:?!. — 7, (ag) 

^ \ ^ >, 0(^,J)•••0(/,A;)0(m,r^)' ^ ^ 

Am 

< prev. > 
[ constr. J 

where we have used the fact that if Xm = A„, then 0(m, n) = f\m, n) + f*{m, n). 

2. Previous constraints do not determine that Xm = A„. Then the situation is as it was 
for the pair (i, j), and the sum will split in the following way 
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E 

\ constr. J 



r(i,...,n) 



f{m,n) 



r(l,...,n) 



f Am.: ■■■ 



<t>{.h3)---<t>{Kk)(l){m,n) 



+nm,n) y: 



r(l,...,n) 



(A7) 



} 



(l>{hj)---(t>{l,k)(t){m,ny 



where in the second term on the right hand side no new constraint has been introduced. 

The geometrical mapping follows because in case (2), every time a factor appears, 
we have a new constraint forcing the two spins to be parallel. On the other hand, in the 
percolation model, such a factor implies that the two vertices are connected, and therefore 
belong to the same cluster. When a factor /* appears, no new constraint is added, so that 
the two vertices can be assigned spins at random. In case (1), on the other hand, the two 
spins are already forced to be parallel by virtue of some previous constraints. According 
to the foregoing argument, this means that they are both connected (at least indirectly) to 
some other spin and therefore that they already belong to the same cluster. In this case it 
does not matter anymore whether these vertices are also linked directly (a factor /^) or not 
(a factor /*). Hence, we see that vertices which belong to the same cluster in the percolation 
picture must be assigned parallel spins while the clusters themselves have randomly assigned 
spin values. 

When all pairs have been covered, the set of constraints of a particular sum specifies 
exactly which particles belongs to which clusters in the original bond-percolation model 
configuration. Since the expression for X^a^ r(l, . . . , n) contains sums over all possible con- 
straints, it can be rewritten as a sum over all possible clusterings of the original percolation 
configuration. Thus, let us define 



P( conn.) 




(A8) 



Then we can write that 
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^r(l,...,n)= Y: E P(conn.) (A9) 

^rn. ( all possible ] f a™, consistent 

{ connectivity \ 1 with the 

I states J 1 connectivity 
[ state 

where the sum over all spins is consistent with the connectivity state in the sense of the 
geometrical mapping, i.e., that all vertices within a single cluster must be assigned the same 
spin. 

If we now assume that two vertices are connected with an unnormalized probability 
and disconnected with an unnormalized probability /* (see section ^V]) , we have the 
interpretation that P(conn.) is the unnormalized probability of finding a specific perco- 
lation configuration of connected and disconnected vertices. Therefore, given a function 
. . . , ra) defined in the percolation model, of the coordinates of the vertices, we can 
define its "average" as 

(F(l,...,n))p= P(conn.)F(l,...,n). (AlO) 

r conn. \ 

\ states J 

Similarly, we can have averages performed in the Potts fluid, which will be denoted by ( ) ^. 
Given a quantity G defined in the Potts fluid, we have 

(^(1, Ai; . . . ; n, A„)), = E ^"(1, Ai; . . . ; n, A„) [] j)- (AH) 

{Am} i>j 



Note that Ai and A2 are assumed fixed. Eq. (|A9|) now implies that in general, for any 
quantity G(l, Ai; . . . ; n, A„) defined in the Potts fluid system, we have that 

(G'(l,Ai;...;n,AO), = ( E G(l,Ai;...;n,A„)\ , (A12) 

HA™: cl} 

where X)x\ ■ n\\ means a summation over all spin configurations which are consistent with 
a given clustering in the sense of the geometrical mapping (i.e., such that all the spins in a 
single cluster are parallel). 

This fundamental relation allows us to relate to each other quantities in the percolation 
model and quantities in the spin model. In particular, let us consider the function fi(i,j) 
defined in the percolation model as 
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1 if i, j belone to the same cluster 

nit,j) = { . (A13) 

otherwise 

If we now denote Pd(l, 2) the probabihty that the two white circles labeled 1 and 2 belong 
to the same cluster, then 

P,Kl,2) = Ml,2))p. (A14) 

We shall now prove that 

{Q{1, 2))^ = lim /-l^ 5: ^{XMX2)) , (A15) 

\ Ai,A2 / g 

where ip{X) equals (s — 1) if A = 1 and —1 otherwise [see Eq. (|1.3|)]. 



By using the fundamental relation Eq. (|A12D , we have 



(j^ E ^(Ai)^(A2)) = E E ^(Ai)^(A2)) (A16) 

\ Al,A2 I g \ {Xm-d.} Xl,X2 ' p 

We now separate the average on the right hand side into two parts: when Ai and A2 belong 
to the same cluster and when they do not. Then, from the definition of the percolation 
average, Eq. (|A10|) , 



( E EV^(AiMA2)) = E ^conn.) E E V^(Ai)^(A2)^^(l,2) 

\{A™.:d.} Ai,A2 / r conn. \ {Am.: ci.} Ai,A2 

\ states J 

+ nconn.) E ^(Ai)^(A2) [1-^1,2)]. (A17) 

r conn. \ {A„:ci.} Ai,A2 

\ states J 

The first sum contributes only if Ai,A2 belong to the same cluster, while the second con- 
tributes only if they belong to separate clusters. Because of the geometrical mapping, 
iplXi) = 'ip{X2) in the first sum. Also, in the sum over {A^ : c/}, every cluster other than 
the one containing Ai and A2 contributes a factor s, the number of possible spin assign- 
ments. The cluster containing Ai and A2 , on the other hand, contributes a factor (s — 1)^ 
if Ai = A2 = 1 and (—1)^ for the (s — 1) other possible choices for Ai = A2. Therefore, 

E E ^{XimWl, 2) = s^-' [is - If + is- n{l, 2) (A18) 

{A„i: d.} Ai,A2 

31 



where K is the total number of clusters in the configuration. Hence, 

E ^^^^^ E E^(^i)^(Wl,2)= E P(conn.).^-^(.-l + lMl,2) 

/ conn. \ ^ {\m.ci.}\iM / conn. \ 

\ states; \ states/ 

= (s^fi(Ai,A2))^. (A19) 

The second term on the right hand side of Eq. ( [A 17] ) contributes only if Ai,A2 belong 
to different clusters. The cluster containing Ai contributes a factor (s — 1) if Ai = 1, and 
a factor (—1) in all the other (s — 1) cases. The same holds for the cluster containing A2. 
The K — 2 remaining clusters contribute each a factor s. Since the values of Ai and A2 are 
assigned independently of each other, we have that 

E E ^(Ai)^(A2) [1 - ^^(1, 2)] = s^-' [{s -!) + (.- 1)(-1)]^ [1 - 2)] = 0. 

{\m.:d.} Ai,A2 

(A20) 



Combining Eq. ([AT9| ) and (|A20D , we see that 



lim/-l^ E V^(Ai)^(A2)) = Ml, 2)),. (A21) 

\ Ai,A2 / g 

Now, separating the case Ai = A2 = o" from the case Ai = cr 7^ A2 = ^, we have that 

(E ^(AiMA2)\ = Y.^\a)R{a,a) + Y.^{a)mR{<y.i). (A22) 

\Ai,A2 I g o- 

where 

i?(Ai,A2)= E r(Ai,A2,A3,...,n), (A23) 

A3,---,An 

in which we have written explicitly the spin variables in the function r(l, 2, . . . , n) for clarity. 
It is important to see that because of the summation over A3, ... , A„, -R(Ai, A2) is independent 
of the specific value of Ai and A2, and depends only on whether they are equal. Because of 
this, we now have that 

Y,i)'^{a)R{a,a) = {s - lfR{a = l,a = 1) + (s - l){-lYR{a = a,a = a), (A24) 
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where a is some arbitrary value of the spin different from 1. The factor (s — 1) in the last 
term on the right hand side represents the possible choices of this spin a 7^ 1. As a result, 



1 



lim cr) = lim -R(a, a) 



(A25) 



where a 7^ 1, but is otherwise arbitrary. 
Similarly, 

i) = {s - Xf a) + 1)] + (s - - 2) r^), (A26) 
where ex ^ Tj are arbitrary values of the spin which are both different from 1. As a result. 



hm ^(a)V^(O R{a, i) = - lim R{a, r^). 

s - 1 s^l 



(A27) 



Combining Eqs. ([A25D and ( |A27|) we have that 

lim/^ Y ^(Ai)^(A2)) =\im[R{a,a)-R{a,r^)] 



(A28) 



Ai,A2 



where a^r] ^ 1 and rj. Comparing with Eqs. ([A14|) , ( [A21| ) and ([A23|) , we finally obtain 
that 



P,^(l,2) = lim 



Y r(a,a,A3, ...,A„) - ^ r(a, 77, A3, . . . , A, 

A3,...,A„ A3,...,A„ 



(A29) 



For the specific example of the four- vertex graph mentioned at the beginning of this section, 
for example, this means that 



ciK 





• • 




CDF 




= lim 




6 6 





o — o o — o 



(A30) 



a a a rj 

This clearly holds for all the diagrams in the CDF expansion. If we restore the integration 
on the black circles, we find that this proves the equivalence, diagram by diagram, of the 
CDF algorithm and the expansion we derived in the present paper. 
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TABLES 

TABLE I. The integrals contributed by triangular diagrams. The first column is the notation 

and the last is the value. This value is the D th power of the one dimensional integral whose 
integrand occupies the second column. The variables x and y are integrated upon. 

h Fi{x) Fi{x - y) Fi{y) (Sa^)^ 

h F2 (x) F2 {x - y) Fi {y) {Ada - a^)^ 

h F2{x) F2{x - y) F2{y) {3d^f 



TABLE IL The integrals contributed by square diagrams with four lines. The first column is 
the notation and the last is the value. This value is the D-th power of the one dimensional integral 
whose integrand occupies the second column. The variables x, y, and z are integrated upon. 

Ji Fi (x) Fi {x - y) Fi (z) Fi {z - y) (f a^) 

J2 Fi (x) Fi {x - y) F2 {z) F2 {z - y) (sda^ - \a^Y 



h Fi{x) F2{x - y) F2{z) F2{z - y) (Od^a - |a3 

■h Fo (x) F2 {x - y) Fhiz) F2(z- y) f if ri'^ 



D 

D 
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TABLE III. The integrals contributed by square diagrams with five fines. The first column is 
the notation and the last is the value. This value is the D-th power of the one dimensional integral 
whose integrand occupies the second column. The variables x, y, and z are integrated upon. 





Fi{y)Fi{x) 


Fi{x 


-y)Fi{z) 


Fi{z 


-y) 




K2 


Fi{y)Fi{x) 


Fi{x 


-y)F2{z) 


F2{Z 


-y) 




K3 


Fi{y)F2{x) 


F2{X 


-y)F2{z) 


F2{Z 


-y) 


(8(Pa - Ada'' + la^y 


K4 


F2{y)Fi{x) 


F2{X 


-y)Fi{z) 


F2{Z 


-y) 


(8da2 - fa^)" 


K5 


F2{y)Fi{x) 


F2{X 


-y)F2{z) 


F2{Z 


-y) 


(ed^a-da' - ^a^)^ 


Kg 


F2{y)F2{x) 


F2{X 


-y)F2{z) 


F2{Z 


-y) 





TABLE IV. The integrals contributed by square diagrams with six lines. The first column is 
the notation and the last is the value. This value is the D-th power of the one dimensional integral 
whose integrand occupies the second column. The variables x, y, and z are integrated upon. 





Fi{x)Fi{y)Fi{z)Fi{x 


-y)Fi{x- 


z)Fi{z 


-y) 




L2 


F2ix) F2{y) F2{z) Fi{x 


-y)Fi{x- 


z)Fi{z 


-y) 


{6da'-2a^f 


L3 


Fi{x) F2{y) F2{z) F2{x 


-y)F2{x- 


z)Fi{z 


-y) 


{8da''-4a^f 


La 


F2{x)F2{y)F2{z)F2{x 


-y)F2{x- 


z)Fi{z 


-y) 


{6d^a - 2da'f 




F2ix) F2{y) F2{z) F2{x 


-y)F2{x- 


z)F2{z 


-y) 


{AdY 
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o 



o 



o 



o 



o 



o oo 



o o oo 







0.2 



0.4 



7] 



0.6 



0.8 



FIG. 1. The percolation threshold, Be, as a function of the aspect ratio for hard core squares 
with a soft shell, in dimension 2. The line represents the theoretical calculation, obtained from 
the AFC algorithm, biased with a value of 7 = 43/18 [19]. The circles are the results of computer 
simulations. 
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FIG. 2. The percolation threshold, Be, as a function of the aspect ratio for hard core squares 
with a soft shell, in dimension 2. The line represents the theoretical calculation, obtained from the 
AFC algorithm, biased with a recently obtained value of 7 = 1.94 [6]. The circles are the results 
of computer simulations. 
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FIG. 3. The percolation threshold, Be, as a function of the aspect ratio for hard core cubes 
with a soft shell, in dimension 3. The line represents the theoretical calculation, obtained from 
the AFC algorithm, biased with a value of 7 = 1.74 [20]. The circles are the results of computer 
simulations. 
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FIG. 4. The percolation threshold, Be, as a function of the aspect ratio for hard core hypercubes 
with a soft shell in dimension 4. The line represents the theoretical calculation, obtained from the 
AFC algorithm, biased with a value of 7 = 1.44 [19]. The circles are the results of computer 
simulations. 
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FIG. 5. The percolation threshold, Be, as a function of the aspect ratio for hard core hypercubes 
with a soft shell in dimension 5. The line represents the theoretical calculation, obtained from the 
AFC algorithm, biased with a value of 7 = 1.2 [19]. The circles are the results of computer 
simulations. 
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